Integrative systems biology analysis of barley transcriptome ─ hormonal signaling against biotic stress

Biotic stresses are pests and pathogens that cause a variety of crop diseases and damages. In response to these agents, crops trigger specific defense signal transduction pathways in which hormones play a central role. To recognize hormonal signaling, we integrated barley transcriptome datasets related to hormonal treatments and biotic stresses. In the meta-analysis of each dataset, 308 hormonal and 1232 biotic DEGs were identified respectively. According to the results, 24 biotic TFs belonging to 15 conserved families and 6 hormonal TFs belonging to 6 conserved families were identified, with the NF-YC, GNAT, and WHIRLY families being the most prevalent. Additionally, gene enrichment and pathway analyses revealed that over-represented cis-acting elements were recognized in response to pathogens and hormones. Based on the co-expression analysis, 6 biotic and 7 hormonal modules were uncovered. Finally, the hub genes of PKT3, PR1, SSI2, LOX2, OPR3, and AOS were candidates for further study in JA- or SA-mediated plant defense. The qPCR confirmed that the expression of these genes was induced from 3 to 6 h following exposure to 100 μM MeJA, with peak expression occurring between 12 h and 24 h and decreasing after 48 h. Overexpression of PR1 was one of the first steps toward SAR. As well as regulating SAR, NPR1 has also been shown to be involved in the activation of ISR by the SSI2. LOX2 catalyzes the first step of JA biosynthesis, PKT3 plays an important role in wound-activated responses, and OPR3 and AOS are involved in JA biosynthesis. In addition, many unknown genes were introduced that can be used by crop biotechnologists to accelerate barley genetic engineering.


Introduction
Cereals are directly related to the food source and energy supply and are always fronted with a wide range of biotic stresses that negatively affect their survival and growth [1,2]. These stresses can cause by living organisms such as bacteria, viruses, fungi, nematodes, and pests [1,3]. Among biotic stresses, fungal pathogens play a fundamental role in causing 70 to 80% of cereal diseases [4,5]. Barley (Hordeum vulgar) is very susceptible to fungal pathogens [6] like Fusarium graminearum and Puccinia hordei, which are the serious threats in wet or semi-wet regions [7,8]. On the other hand, pests are other destructive pathogens of cereals, which reduce plant growth without specific leaf symptoms [9]. Hence, identification of the key genes related to the resistance is the first step in development of new breeding lines to manage the diseases [10][11][12]. For example, the breeding for barley resistance to Rhopalosiphum padi has revealed that several resistance genes are associated with the decrease in aphid growth [13]. Several gain-of-function and loss-of-function mutants have shown the genes related to fundamental defense and effector-triggered immunity [14]. Today, large-scale data generated from functional genomic studies on plant-pathogen interactions have provided massive information [15]. Therefore, access to OMICS information on the behavior of defense genes against pathogens helps to decipher the gene networks that led to the identification of pathogen specific-responses [16]. In recent years, computational systems biology has collected exquisite opportunities to dominate biological complexity [17,18]. Systems biology approaches, including meta-analysis and gene network analysis, are wellknown strategies for integrating these data and discovering particular molecular interactions [19,20]. These approaches provide a robust statistical framework for re-evaluating key findings, improving sensitivity by increasing sample size, testing new hypotheses, and characterizing novel gene candidates for plant breeding programs [21][22][23][24].
Therefore, progress in whole-genome transcriptome analysis allows the analysis of gene regulation and plant-pathogen interactions in various conditions, allowing the identification of pathogen-specific biomolecular networks and hormonal signaling pathways. For example, the important responses to pathogens may be due to the cross-talks among hormone-related pathways [25] such as ethylene (ET), salicylic acid (SA), jasmonic acid (JA), abscisic acid (ABA), and/or signaling due to calcium [26], activated oxygen species (ROS) [27], and phosphorylation chains [28]. Although many signaling pathways interact, regulatory factors such as transcription factors (TFs), microRNAs (miRNAs), and protein kinases (PKs) regulate plant defense responses. Most plant immune responses are transcriptionally controlled [29,30] and some pathogen-derived factors disrupt host transcriptional complexes to increase their virulence, highlighting the importance of transcriptional regulation of host defenses [31].
In defense responses, TFs play distinct roles in modulating gene expression and activating signaling cascades [32,33]. It has shown that the WRKY family, which is considered the core of the plant immune system, plays an essential role in phytohormone signaling and pathogen defense [34]. Members of MYB/bHLH regulate distinct cellular processes, such as responses to biotic stresses, hormonal signaling, and hypersensitive reaction (HR). In addition, the bZIP family plays an important role in regulating pathogen-responsive pathways and mitigating pathogens [35].
Plants are equipped with an array of defense strategies against pathogens [36]. Systemic acquired resistance (SAR) and induced systemic resistance (ISR) are two forms of induced resistance pathways [37]. ISR is a nonspecific resistance mechanism that is activated by infection, and acts against a considerably broader spectrum of both (hemi-) biotrophic and necrotrophic pathogens [37,38]. On the other hand, SAR in systemic and uninfected tissues provides long-term protection against (hemi-) biotrophic pathogens. In addition, SAR is mostly studied as a leaf-to-leaf response [39] and is activated during the local defense response [39], while ISR is triggered by the colonization of roots by certain non-pathogenic rhizosphere microorganisms such as plant growth-promoting rhizobacteria (PGPR) and plant growth-promoting fungi [37][38][39][40].
Many hormonal signaling are involved in both defense systems. For example, SAR requires the synthesis of SA, which triggers the expression of a well-known set of pathogen-related (PRs) genes [41,42], while ISR is dependent on JA and ET signaling pathways [43]. These accumulating signaling molecules moderate the defense responses, and when they are used exogenously, they are enough to induce the resistance [39]. However, ISR and SAR together present a better protection than each of them alone, and this shows that they can act increasingly in inducing resistance to the pathogens [39]. Furthermore, plants rarely encounter only one stressor in a natural habitat when they evaluate defense responses. Therefore, understanding the interactions between SAR/ISR and other signaling cascades help to raise our knowledge for strengthening the plant defense against biotic stresses.
Many studies on SA-mediated plant defense have demonstrated that SA is not only a major regulator of SAR and induces the PRs, but also this hormone is necessary to induce the HR [44]. There is a complex network of SA signaling mediated by the NPR1 and mitogen-activated protein kinases (MAPKs) [45]. The HR control by TFs such as MYB depends on SA accumulation, but not on NPR1, a signaling protein downstream of SA accumulation [46]. Moreover, SA-dependent defense responses are effective against biotrophic pathogens that auto-activate the HR [46]. Further, it has shown that the suppression of SSI2 reduces the 18:1 fatty acids and induces phenylalanine ammonia-lyase-mediated SA accumulation, which auto-activates the HR [47].
To gain a clear insight into the mechanisms involved in the response of barley to various pathogens and hormonal signaling, we used meta-analysis and co-expression gene network analysis of barley transcriptome. We identified some key genes and significant gene networks in response to pathogens and hormones. These genes were identified through gene enrichment analysis of metabolite pathways, TFs, PKs, and microRNAs.
For instance, TFs interact with enhancers to coordinate gene expression transcriptionally [48], while miRNAs are gene regulators transcriptionally or post-transcriptionally [49]. Since both regulators demonstrate great impact on plant genetic systems, the circuiting of miRNAs-TFs will allow the orchestration of diverse biological processes with high reliability [49,50]. Therefore, observing and understanding a dynamic relationship among these regulators is interesting. In addition, some significant hormone-and pathogen-responsive hub genes were validated after methyl jasmonic acid (MeJA) treatment. Plant protection could be facilitated by the exogenous application of the MeJA, which has demonstrated great potential as a defense inducer for fungi and pest attacks [51]. The MeJA stimulates induced systemic SA/JA/ET pathways in the infected tissues, thus inducing the expression of genes related to the resistance. MeJA affects many genes related to both SAR and ISR pathways. Previous studies have shown that PR1, AOS, and LOX2 respond positively to the MeJA treatment, protecting the plant from infection [52,53]. Therefore, more accurate information on barley-pathogen interactions might be helpful for plant biotechnologists to supply the pathogen-tolerant cultivars.

Materials and methods
We integrated the transcriptome datasets of barley through systems biology approaches to identify responsive key genes involved in the pathogen-hormone interactions. First, the metaanalysis discovered the DEGs, and then the co-expression gene network analysis grouped the DEGs into the modules (Fig 1).

Data collection
The raw microarray expression data of barley exposed to diverse biotic stresses, including fungi, bacteria, viruses, and insects, or various hormonal treatments, including JA, gibberellic acid (GA), and MeJA, were retrieved from publicly available databases of Gene Expression Omnibus (GEO), Sequence Read Archive (SRA), and Array Express. The data selection was done based on the Minimal Information about a Microarray Experiment (MIAME) requirements [54]. Therefore, the studies with relatively similar genetic backgrounds without any mutated or transgenic samples were chosen with suitable quality and biological and technical replicates of controls and treatments [55]. The total data included ten studies with 479 biotic stress samples and four studies with 46 hormonal treatment samples (S1 Table). The raw data included two platforms, Affymetrix (Accession: GPL1340) and Agilent (Accessions: GPL14904 and GPL15513). The probe-gene maps and probe-annotation files were retrieved from the Affymetrix database.

Pre-processing and normalization
The normalization of the Affymetrix data for each array was performed with Robust Multichip Average (RMA) algorithm [56] by the Expression Console software (V. 1.3.1). Briefly, after background correction of this platform by the RMA method, the Quantile method normalization was carried out, and then the Median Polish method was applied to summarize the probe sets [57]. Two color Agilent data were processed using Flex Array software (V. 1.6.3). Briefly, the background correction was performed using the Normexp method to estimate the intensity of expression values [58], and then the normalization was performed using the Quantile method [59]. For the normalization of one color Agilent data, we applied the LIMMA package and LOESS method in the R program [60].

Meta-analysis and identification of DEGs
At the first, the control and treatment samples of each study were separately defined. Then the guardianship data, including TAIR IDs, was used for meta-analysis using the MetaDE R package. The merged datasets were filtered again except for 20% of the expressionless genes (with low expression intensity) and 20% of the non-educational genes (with quantitative changes).

PLOS ONE
Finally, the RankProd method, a non-parametric method, was applied to genes rank based on the fold change (FC). Principally, the RankProd algorithm calculates pairwise FC with replicates for each gene between treatment and control samples in both directions, respectively, and converts FC into ranks among all genes, then searches for genes that are continuously topranked across replicates [61]. Converting FC into ranks overcomes the heterogeneity among multiple datasets and an overall ranked gene list is produced based on the P-value and the False Discovery Rate (FDR) of each gene. In this method, up-and down-regulated genes with log2 FC >1 or log2 FC <−1 and FDR � 0.05 were considered significant differentiallyexpressed genes (DEGs) [62].

Gene enrichment analyses
Gene ontology (GO) enrichment analysis of the DEGs to finding of the molecular pathways responsive to pathogens and hormones was performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) and DAVID resources [63,64]. GO was performed on the unique list of TAIR IDs of the DEGs. The results were categorized in Molecular Function (MF), Biological Process (BP), and Cellular Component (CC). GO terms with the corrected P-value � 0.05 were considered significant [65].

Identification of TF, PK and miRNA families
The DEGs were BLASTed against the iTAK database [66] to identify the TFs and PKs. Next, the psRNATarget database [67] was used to identify candidate microRNAs. Then, the genes targeted by miRNAs were considered to GO analysis. The WEGO software [68] was used for plotting GO results, and the Agrigo software was used to enrich the GO terms.

Gene co-expression network analysis and identification of modules
In order to identify the DEGs with similar expression pattern, a Weighted Gene Co-expression Network Analysis (WGCNA) in R space was performed on the matrix of normalized expression values. For this purpose, the meta-analysis output file (including DEGs) was merged with the initial normalized data and then the treatment and control samples were separated in the columns. First, a similarity matrix was calculated based on Pearson correlations between each DEG pair and converted into an adjacency matrix by applying the soft-thresholding β power. This power makes the adjacency matrix to be the continuous value between 0 and 1 and the co-expression similarity was raising to achieve scale-free topology. Following, the Topological Overlap Matrix (TOM) [69, 70] was calculated for hierarchical clustering analysis. Then, to discover modules, a dissimilarity matrix is obtained (dissTOM) that resulted in genes with similar expression clustered in the same gene module. Finally, a dynamic tree cut algorithm was used to identify co-expression gene modules. This method can detect nested clusters and have been shown to better detect outliers [71]. The major parameters were defined with a minimum size cut-off of 30 genes and a threshold of 0.3 and a minimum size cut-off of 100 genes and a threshold of 0.99 for avoiding abnormal modules in dendrogram, on the biotic and hormonal data respectively. As a result, an undirected weighted network with scale-free topology composed of modules of barley DEGs with correlated expression during pathogen infections was obtained.

Identification of hub genes
Network visualization and calculation of topological properties for each module were performed using the Cytoscape software (V. 3.6.1). To identify the hub genes, computational algorithms of Maximal Clique Centrality (MCC) were applied using a plug-in of Cytoscape, CytoHubba [72, 73] followed by GO enrichment analysis by DAVID database.

Promoter analysis
To find the upstream flanking region of DEGs (1000 bp) [74], a BLAST against the Ensemble plant database was performed [75]. Then, the MEME [76] database was used to identify the conserved cis-acting elements. The motifs were counted, and the abundance of regulatory elements in the promoter region of various genes was investigated. Next, we used Tomtom v 5.0.1 tool [77] to remove redundancy motifs and to determine known CRE based on the motif database of JASPAR (https://jaspar.genereg.net/) [78] with a threshold E-value of 0.05. GOMO tool [79] was also used to identify the biological roles of the cis-acting elements.

Protein-protein interactions
The protein-protein interaction network (PPI) was constructed based on the both datasets (biotic stress/hormonal treatment). The unique list of TAIR IDs of each dataset was BLASTed against the STRING database with default parameters (with the highest required interaction score = 0.9). The highest confidence was selected to make more comprehensive network and to study more important connectivity's. Finally, the network was drawn by Cytoscape software and key genes in networks were screened out by the score. In this network, high degree could represent essential genes [80]. In the network, nodes represent genes and edges represent the interactions between the nodes. The highly important nodes in the network, the hub nodes, were obtained according to the k-core.

Plant materials and growth conditions
The seeds of hordeum vulgare genotype Oxin were sown in trays, germinated, and grown in a glasshouse in the pots that contained peat moss, perlite vermiculite, and soil (1:1:3) under a 16-h light:8-h dark cycle at 25˚C (S1 Fig).

Stress conditions and sampling
The seedlings in the 3 leaf-stage were sprayed and irrigated with 100 μM MeJA plus 0.1% Tween-20 in 0.1% ethanol. For the controls, 0.1% Tween-20 solution in 0.1% ethanol was used. Plastic covered all the pots. Leaf samples were collected at 3, 6, 24, and 48 hours after the treatment in three biological replicates. For RNA extraction, the leaves were immediately frozen in liquid nitrogen and stored at −80˚C.

RNA extraction, DNase treatment, and cDNA synthesis
Total RNA was extracted from leaf samples using a Column RNA Isolation Kit (DENAzist Asia Co., Mashhad Iran) according to the manufacturer's instructions. The quantity and concentration of RNA were measured using a NanoDrop ND 1000 Spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). The integrity and quality of RNA were checked by visual observation of 28S and 18S rRNA bands on a 1.2% agarose gel. Before use, RNA samples were stored at −80˚C. DNase treatment of RNA was carried out using the RNase-free DNase kit (Thermo Fisher Scientific, Waltham, Massachusetts, USA) according to the manufacturer's instructions. The treated RNA was rechecked by NanoDrop and agarose gel. Then, 1 μg of DNase-treated RNA was used for first-strand cDNA synthesis using a SinaClon BioScience kit (Karaj, Iran) according to the manufacturer's instructions. The cDNA samples were stored at −20˚C prior to use.

Selection of candidate genes for gene expression analysis
To validate the transcriptome data analysis, 6 genes with the highest interaction in the PPI network related to hormone signaling and disease resistance-related pathways were selected. These candidate genes are the most important hub genes that show interactions between SAR and ISR. These genes play a central role in defense responses and hormone signaling [81]. The first candidate gene is SSI2 (AT2G43710), a hub gene in the biotic turquoise and brown modules, and hormonal blue module, which interacts with fatty acids (FAs) that affect defense signaling (S5 Fig). The SSI2 represses plant defense by affecting the generation of a lipid molecule necessary for SA synthesis and the NPR1-independent pathway. Moreover, the SA and JA signaling-related regulatory genes inducing the PRs are activated in response to numerous pathogens and ultimately enhance resistance [82]. Among the PRs, PR1 has been frequently used as a marker gene for SAR in many species [83]. As the second candidate gene, we selected PR1 from the biotic brown module. The third is PKT3 (AT2G33150) which is involved in longchain fatty-acid beta-oxidation prior to gluconeogenesis during seed germination and seedling growth [84]. This enzyme was identified in the biotic brown and hormonal green modules that confers sensitivity to 2,4-dichlorophenoxybutiric acid (2,4-DB). Required on local and systemic induction of JA biosynthesis after wounding and pathogen infection and involved in JA biosynthesis during senescence [85].
Other candidate's key genes, such as AOS, OPR3, and LOX2, showed a high positive correlation and interaction with each other. These genes were identified in biotic blue or hormonal turquoise modules, hormonal red modules, and biotic turquoise modules, respectively (S5 Fig). The induction of the JA-responsive marker gene by different JA derivatives was equally sensitive to SA-mediated suppression. Activation of genes encoding key enzymes in the JA biosynthesis pathway, such as LOX2, AOS, AOC2, and OPR3, was also repressed by SA, suggesting that the JA biosynthesis pathway may be a target for SA-mediated antagonism [81]. These show that JA biosynthesis is under the control of a positive feedback regulatory system and that the down-regulation of JA biosynthesis-responsive genes by SA may be a critical mechanism in SA-JA cross-talks [81].

Primer design
Primers were designed using Allele ID 7 and Vector NTI 11 software for the reference and candidate genes ( Table 1).

Validation of candidate genes
First, primer specificity was confirmed by PCR and sequence analysis. To minimize pipetting error, the cDNA samples were diluted 1:20 by using nuclease-free water, and 5 μL cDNA was used for qPCR. Relative qPCR was performed in a 20 μL volume containing 5 μL cDNA (diluted), 10 μL RealQ Plus 2x Master Mix Green (Sinuhebiotech, Shiraz, Iran), and 0.7 μl of primers (10 μM). The amplification reactions were carried out in a Line-gene K thermal cycler (Bioer, China) under the following conditions: 15 sat 95˚C, 45 cycles of 94˚C for 15 s, Ta temperature for 15 s, and 72˚C for 20 s. After 45 cycles, the specificity of the amplifications was tested by melting curve analysis by heating from 50 to 95˚C. All amplification reactions were repeated three times under identical conditions and included a negative control. In addition, 6 serial dilution samples (1:10) provided using specific PCR products of a cDNA template for each primer were used in qPCR as positive controls. Besides, we used these standards to calculate the qPCR efficiency.

qPCR data analysis
The relative expression was calculated according to the 2 ─ΔΔCT method [86]. The 6 C T values for each sample were calculated using the Line-gene K software. Three technical replicates were run for each sample in qPCR. The actin gene (AT3G18780) was used as the reference gene for data normalization [87,88].

Statistical analysis
Analysis of variance followed by Duncan's multiple range test was performed using MINITAB (Minitab, Inc., Pennsylvania, USA). In all cases, differences were regarded to be statistically significant at a P-value � 0.05. All experiments were performed in triplicate, illustrated using the GraphPad prism software (GraphPad Software, Inc., San Diego, CA).

Pre-processing and normalization
Pre-processing was applied to the raw data to reduce the impact of existing noises. Since reducing heterogeneity among studies is a necessary step for the integration of the data, the normalization methods of RMA and LOESS were applied to Affymetrix and Agilent platforms, respectively. The proposed framework of the normalization reduced the batch effects and facilitated the direct merging (S2 Fig).

Meta-analysis and identification of DEGs
In this study, the RankProd method was used to identify DEGs for biotic stress and hormonal treatment studies (S1 Table). The results of the meta-analysis may contain duplicated genes occurring by multiple probes mapped into the same gene symbols. These duplicated genes, if any, are filtered out by different methods after obtaining gene values for each dataset. For each pair of duplicated genes, the up or down-regulation was considered based on the smaller Pvalue. Thus, we found that from 44 duplicated genes in the biotic data, 10 genes remained in the form of up-regulated and 34 genes remained in the form of down-regulated, as well as from 4 duplicate genes in the hormonal data, 2 genes were considered as up-regulated and 2 genes considered as down-regulated.
The KEGG pathway analysis for biotic DEGs revealed significant enrichment in pathways associated with the biosynthesis of secondary metabolites, antibiotics, and amino acids ( Fig 3B).
In addition, the most percentage of hormonal DEGs ( Fig 4A) were involved in the oxidation-reduction process (15% genes), defense response (5.5% genes), and the responses to oxidative stress (3.2% genes) in BP. In MF, metal ion binding and oxidoreductase activates were dominant. The category of CC and KEGG pathway analysis was the same as the components of biotic DEGs ( Fig 4B).

Co-expression analysis and module identification
The R package WGCNA was applied to 1232 biotic and 308 hormonal DEGs to identify coexpressed modules. The soft-threshold power was one of the main parameters in WGCNA, making the constructed network more in line with the scale-free network characteristics and average connection degree of the gene co-expression module. When the soft-thresholding β power was 10, the scale-free topology fitting index (R2) was optimized to � 0.8, and this power was selected to build a hierarchical clustering tree. By the dynamic tree-cutting algorithm, DEGs were grouped into 6 modules for biotic data ranging from 81 to 431 genes per module (S3A Fig and Table 2) and 7 modules for hormonal data ranging from 141 to 926 genes per   Fig and Table 3). Finally, each module was assigned on a unique color label, which was used as an identifier for the analyses. The hierarchical clustering was illustrated in the following. As the dendrogram demonstrated, we achieved two clusters, including 4 sub-clusters for biotic data ( Fig 5A) and 5 subclusters for hormonal data ( Fig 5B). In most modules, including blue, brown, and yellow or red and green modules, biotic DEGs were expressed similarly according to Multidimensional Scaling (MDS) (Fig 5C). In addition, hormonal DEGs in the modules of blue, green, and turquoise or brown, red, and yellow exhibited a similar expression pattern ( Fig 5D). The brown and blue modules showed the strongest gene-gene interconnected based on TOM dissimilarity distances.
In the eigengene adjacency heatmap, the slope of the variance in color from black or blue to yellow represents the connectedness of genes for different modules from strong to weak, whilst the red highlights are the significant modules related to biotic and hormonal data (Fig 5A and  5B). The brown and blue modules exhibit the strongest gene-gene interconnected based on TOM dissimilarity distances. The adjacency heatmap of eigengenes module represents the summary of the results in the form of color-coded tables, illustrating the good agreement between consensus modules of each dataset or between specific and common modules of biotic and hormonal datasets. It reveals that most modules, especially the brown and blue modules, are preserved in biotic and hormonal treatment. Moreover, based on the TOM, light-yellow color represents a higher overlap between modules, and a progressively dark color represents a lower overlap between modules. Blocks of light colors along the diagonal are the modules (Fig 5E and 5F).
Finally, the gene enrichment analysis was performed for BP category to understanding the biological functions associated with the modules (S3 Table). Based on the BP results of biotic modules, the turquoise module was significantly enriched in oxidation-reduction processes, and the blue module was enriched in intracellular protein transport and photosynthesis. The brown, yellow, green, and red modules were related to cadmium ion response, translation, protein ubiquitination, and stress response, respectively.
According to the GO results of hormonal modules (S3 Table), the most important BP terms enriched in turquoise and blue modules were significantly enriched in photosynthesis, response to cytokinin, protein transport, response to cadmium ion, protein folding, isopentenyl diphosphate biosynthetic process, oxidation-reduction process, and ribosome biogenesis. The brown module was related to the defense response to a bacterium, reactive oxygen species metabolic process, oxidation-reduction process, response to wounding, and metabolic process. In addition, 600 common genes (17.1%) were identified between biotic and hormonal modules in which the turquoise (195 genes), blue (169 genes), and brown (117 genes) modules involved the most common genes.

Hub genes
The top 30 hub genes were chosen for each module. GO analysis was performed for 177 biotic and 198 hormonal hub genes. The hub genes related to biotic modules were strongly enriched in photosynthesis, translation, and ribosome biogenesis pathways ( Fig 6A).
Among the hormonal modules, the top enriched BP terms were associated with the oxidation-reduction process, response to cadmium ion, and ubiquitin-dependent protein catabolic process (Fig 6B). In the following, 10 common hub genes (2.7%), such as NF-YC2, SSI2, ARAC1, and PFD1, were identified between biotic and hormonal data (S4 Fig). Besides, hub genes with unknown functions, including AT3G26670, AT1G67620, AT1G05720, AT1G52600, AT5G59140, and AT3G22290, were identified as potential candidates for further investigation.

Identification of transcription factors
To identify important TFs, DEG sequences were used to the BLASTX against the iTAK database. All 24 TFs found for biotic DEGs belonged to 16 different TF families and were either directly or indirectly involved in hormonal signaling and response to biotic and abiotic stresses ( Fig 7A and S5 Table). Members of the auxin/indole-3-acetic acid (Aux/IAA), GCN5-related N-acetyltransferases family (GNAT), nuclear factor Y, subunit C (NF-Y), and WHIRLY families were the top classes. Among these TFs, only 5 families, including bZIP and GNAT, were significantly up-regulated, while most family members were down-regulated. These TFs were more involved in response to wounding, bacterium, karrikin, oxidation-reduction, ABA, JA, and in the arginine biosynthetic process. We found a higher abundance of TFs in the turquoise module (8 TFs) as compared to the brown (7 TFs), green (4 TFs), blue (2 TFs), and red (2 TFs) modules. The WHIRLY (in the turquoise module), GNAT, and NF-YC families (in the green module) were also considered as the hub genes (S5 Table).
Furthermore, 6 TFs related to hormonal DEGs were identified, belonging to 6 different TF families (Fig 7B and S5 Table). Among hormonal TF families, only one family of OVATE family proteins (OFPs) was significantly up-regulated, while other families were down-regulated. These TF families were more involved in response to oxidative stress, wounding, bacterium, ABA, and JA signaling pathways. We further extended the study to detect hub genes and TFs in each module. The brown module was enriched for Tify and OFP TF families, whereas the turquoise, red, and black modules included at least one member of GNAT and MBF1 families.
Among the identified TF families of biotic or hormonal studies, there was only one common significantly up-regulated TF of TIFY10A (AT1G19180) that was involved in response to JA, defense responses to the bacterium, and regulation of defense responses. relationships, the progressively more saturated blue and red colors indicate the high co-expression interconnectedness. c-d) Multi-dimensional Scaling (MDS) map of each gene located in a module of the c) biotic or (d) hormonal studies. MDS plot demonstrating the similarity of gene expression patterns between different modules. Genes of different modules are marked in different colors. e-f) show network heatmap plots. Co-expression network modules of (e) biotic or (f) hormonal DEGs are presented. In the TOM plot, darker red color represents low overlap as progressively light color represents higher overlap among DEGs. Blocks of darker colors along the diagonal correspond to the modules. Gene dendrograms and module assignments are also displayed on the left side and the top. https://doi.org/10.1371/journal.pone.0281470.g005

Identification of PKs
PKs are important signaling regulators in response to environmental stresses [91]. To identify PKs, DEG sequences were used to the BLASTX against the iTAK database. Three upregulated PKs were identified among the biotic DEGs and classified into RLK (2 types of receptor-like kinase-Pelle (RLK-Pelle), and tyrosine kinase-like kinase (TKL) (families (Table 4). The RLK-Pelles, the largest group, were commonly identified in turquoise and blue modules. In addition, TKL-Pl-4 was identified in the turquoise module.
Among hormonal DEGs, three PKs were identified, which were classified into RLK-Pelle and TKL families (Table 4). Among the identified TKLs, there was one up-regulated (AT5G58950) gene at the red module and the others were down-regulated. RLK-Pelle families were also identified in the black module. There was only one common up-regulated PK detected in the red and turquoise modules among biotic and hormonal PKs. This PK is TKL-Pl-4 (AT5G58950).

Identification of miRNAs
The psRNATarget tool was used for predicting potential miRNAs that target the DEGs. Among the biotic DEGs, 432 miRNAs belonging to 56 conserved families were found (S6 Table). Furthermore, 85 miRNAs belonging to 32 conserved families were identified among hormonal DEGs (S6 Table). It was found that 15 target genes were common across both datasets; of which miR6192 was the most frequent ( Fig 8A).

Cis-acting elements
The 1000 bp upstream flanking region of DEGs was used to predict conserved cis-regulatory elements (CREs). For each biotic and hormonal DEG, MEME analysis identified 11 significant cis-acting elements (Tables 5 and 6). Among the top motifs for biotic and hormonal DEGs, motifs 1 and 4 showed the highest frequency. Moreover, most of the cis-acting elements have related to C2H2 zinc finger motifs like A1267.1, MA1267.1, MA1267.1, and MA1268.1 for the biotic DEGs, as well as MA1267.1 and MA0415.1 for the hormonal DEGs (S7 Table). The GOMO analysis for the motifs found by MEME detected many interesting BFs (Tables 5, 6 and S8). GO term analysis of biotic DEGs indicated that these motifs were involved in the regulation of transcription, DNAdependent, and signaling pathways. WRKY families targeted the main biotic motifs, whereas 85% of the DEGs involved at least one WRKY binding site. Moreover, these motifs were involved in MFs of TF and protein serine/threonine kinase activities (Table 5). GO term analysis of hormonal DEGs showed that the motifs were mainly involved in the regulation of transcription ( Table 6). The key TF families comprised WRKYs, MYBs, bHLHs, and bZIPs. Finally, this analysis identified the motifs related to SA (GO: 0009751) and auxin pathways (GO: 0009733) (S9 Table).

Protein-protein interactions and selection of key genes
It is important to find genes that are regulated for any type of biotic stress or hormonal treatment. We selected 1232 biotic DEGs and 308 hormonal DEGs. The DEGs also represent common responses to defense and may be useful in describing general responses to the biotic stresses in barley and Arabidopsis. For this purpose, the analysis of protein-protein interactions (PPIs) using the STRING database and Cytoscape software was performed (Fig 9A and 9B). The network was plotted for 265 key biotic DEGs, which were involved in the most important mechanisms responsive to the biotic stresses, and 67 hormonal key DEGs, which were related to the hormonal signaling pathways (Table 7).
A list of genes that showed the most protein-protein interactions was prepared (S10 Table). It is interesting to notice that some up-regulated biotic or hormonal hub genes may play a key role in response to fungal pathogens or hormonal signaling. Finally, 4 biotic key genes, including PKT3, which was active in the local and systemic induction of JA biosynthesis after pathogen attack, PR1, which involved in response to a variety of pathogens and a useful molecular marker for the SAR signaling, SSI2, which was active in SA-and JA-mediated defense signaling, and LOX2, which was active in the pathogen resistance and synthesis of JA, were identified. In addition, there were two hormonal key genes, including OPR3 and AOS, which involved in the biosynthesis of JA and other oxylipin signaling molecules, as well as the response to fungal pathogens.
Finally, a PPI network was illustrated to study key genes that are commonly expressed among biotic and hormonal DEGs. We showed a total of 61 common DEGs such as SSI2, PAL, JAZ1, APE1, ribosomal protein L24 (RPL24), and some others significant DEGs especially with unknown function which are simultaneously active in the pathways of resistance to the pathogens and signaling of hormones (Fig 10).
Validation of candidate genes. Here, a significant number of all DEGs were regulated by JA/SA in the same or opposite direction. In this study, 238 co-induced JA/SA DEGs were defined as wide-spectrum pathogen resistance DEGs activated by bacteria, fungi, viruses, and aphids (S10 Table and Figs 9 and 10). To experimentally validate the reliability of the DEGs obtained from the integrative transcriptome analyses, a total of 6 hub genes including SSI2, PR1, PKT3, AOS, OPR3, and LOX2 were candidates for qPCR (Figs 11 and S6). These genes were chosen based on being highly differentiated in response to the MeJA and pathogens.
The candidate genes (P-value � 0.05) responded positively or negatively to the application of MeJA. However, the expression level of SSI2 reached the highest peak at 3 h, about 2 times higher than the control. Notably, SSI2 showed a dramatic suppression within 6 h after treatment, and it showed no significant increase after 24 to 48 h compared to the control, indicating that the SSI2 may perform different functions in response to MeJA signaling (Fig 11A).
Conversely, we found that the expression of PR1 and PKT3 were significantly increased (Fig 11B and 11C). Although the FC at the time point of 48 h slightly decreased compared to the rest of the periods, it showed a significant increase in the expression compared to the control. Meanwhile, for the AOS, OPR3 and LOX2 expression, we found that they were significantly induced in response to MeJA treatment after 24 h (Fig 11D-11F). Additionally, the expression of AOS, OPR3 and LOX2 was induced at 3 to 24 h after treatment, while it decreased to its basal level (restored to normal condition) after 48 h (Fig 11D-11F).

Discussion
In spite of many researches, there has not been sufficient clearance of the molecular mechanisms of defense in barley. For instance, what are the main mechanisms of hormonal signaling against different pathogens or what are the main gene networks? The integrative analysis of the transcriptome is one of the most efficient approaches to finding the answer to biological questions, thus enabling researchers to achieve the huge amounts of valuable information. In this study, based on the biotic stress and hormonal treatment data, we discussed the defense responses of barley and several specialized pathways related to pathogens. In addition, the responsive specific novel genes to defense hormones in barley that are the results of our research were elucidated.

Common and unique biological processes responsible for pathogen resistance and JA signaling
To identify the main DEGs and key pathways, the integration of transcriptomic data retrieved from various biotic stress and hormonal treatment studies in barley, was applied by meta-analysis. Pathogen resistance and hormone signaling are studied using this approach, which allows us to understand the response pathways between pathogens and hormones. The meta-analysis identified 1232 biotic and 308 hormonal DEGs (Fig 2A and 2B and S2 Table). Of 61 common  DEGs between biotic and hormonal studies, 1.2% and 0.8% were categorized as up-and down-regulated, respectively (Fig 2B).
The gene enrichment analysis suggested that DEGs were enriched in various BPs (Figs 3  and 4). The terms of photosynthesis, response to cadmium ion, protein folding, oxidationreduction process, response to CK, defense response to the bacterium, defense response to fungus, cellular response to oxidative stress, response to ABA, JA mediated signaling pathway,

PLOS ONE
response to reactive oxygen species, response to wounding, and response to biotic stresses and hormones, were also found to be significant. For example, PAL is involved in the phenylpropanoid pathway as the primary biosynthetic enzyme and in SA biosynthesis during SAR induction [92,93]. The importance of SA as a key signaling hormone in resistance to diseases has considerable interest. Since SA accumulation is associated with the expression of defense marker genes, the PAL pathway produces the intermediates required for the biosynthesis of not only SA but also many other secondary metabolites, such as lignin [94]. PAL transgenic crops were more resistant to brown planthopper (BPH) [92]. Hence, the expression of some PALs was significantly up-regulated under pathogen infestation, the up-regulated PALs may be correlated with BPH-elicited phenolamide and flavonoid accumulation in crops. Carbohydrate metabolism has been widely shown as an essential pathway influenced by biotic stress responses [95]. Growing organs and tissues are the complex sinks of carbohydrates attracting photosynthesis energy produced by leaves. The source and sink relationships allow the allocation of carbon during a variety of stresses to improve plant performance [96]. The disruption of source-sink connections has been associated with the early pathogenic mechanisms of diseases in plants [97]. The dysregulation of this pathway at the transcriptomic level may be associated with a general plant stress state.
Plant elicitors are attractive and promising factors for plant resistance because they increase host endogenous immunity without direct toxicity to the pathogens [98]. Screening strategies for these elicitors could be mentioned as pathogen-responsive-gene-based (PRGB) approaches. Here, 4% of hormonal (most JA/SA) induced DEGs are commonly regulated by biotrophic and necrotrophic pathogens, which were confirmed by microarray datasets, and were defined as broad-spectrum disease resistance genes (Fig 2). These DEGs are involved in nitrogen storage, photosynthesis, auxin transport, and gibberellin biosynthesis. It suggested that plants under treatment of the elicitors or attack by pathogens have prioritized appropriate defense response over the growth by activating the shared SA and JA responses. The results showed that the JA-responsive marker gene PDF1.2 and SA-responsive marker gene PR1 were uniquely upregulated by JA and SA, respectively. The promoter of PDF1.2 responded to MeJA but not to SA [99]. Thus, certain downstream genes of SA signaling might be free from the impact of JA signaling activation. Therefore, plants simultaneously activate JA/SA signaling pathways to suitable responses to the different pathogens [100]. The allene oxide cyclase (AOC4), JA carboxyl methyltransferase (JMT), AIM1, and PKT3 families are the regulators of seed germination and subsequent seedling growth. Moreover, they have required for the local and systemic induction of JA and MeJA biosynthesis after wounding. It seems to be involved in JA biosynthesis during senescence [101][102][103][104].
The study of the hormonal responses to pathogens is vital because the hormonal signaling pathways regulate defense responses antagonistically or synergistically [105]. It is observed that several genes (such as WRKY families) involved in auxin, ET, JA, and SA signaling were usually affected against all types of pathogens. Fungal pathogens up-regulate ET, benzyl-adenine (BA), and SA signaling pathways while mostly repressing JA-related genes. For example, Erwinia amylovora up-regulates ET and GA signaling genes, while viruses affect some essential genes involved in ET and auxin signaling [105]. Our results confirmed that all barley pathogens deeply affected the hormone-related pathways.

TFs, regulatory factors in response to pathogens
To present the information regarding how genes are regulated in barley against pathogen attacks and hormone treatments, we identified more important TFs. TFs only not act as the activators for the expression of stress-inducible genes but also play a significant role in signal transduction pathways in biological networks and promote plant adaptation to various environmental conditions [106]. Finally, we identified 24 biotic TFs belonging to 25 TF families and 6 hormonal TFs belonging to 6 TF families (S5 Table). The major members of TF families, such as GNAT, NF-YC, and WHIRLY, were the top classes and were also considered hub genes (Fig 7). In this analysis, WHIRLY was up-regulated. Plants have a small family of singlestranded DNA-binding proteins called WHIRLY, which are also involved in defense gene regulation [107]. In barley, WHIRLY binds to the promoter of the HvS40 senescence-related gene, which is induced by pathogens, senescence, and hormones, such as SA, JA, and ABA [107,108]. Transgenic barley with suppressed WHIRLY1 remained almost green phenotype despite having higher SA than the wild-type [109]. The plant-specific WHIRLY is induced by SA and required for both SA-dependent disease resistance and SA-induced SAR gene expression. WHIRLY is required for both full basal and specific pathogens resistance responses. It has been found that NPR1 is required for SAR, and that is an important component of the SA signaling pathway that is activated independently of the key SAR regulator NPR1. However, both NPR1 and WHIRLY1 are required for optimal SA-induced defenses. WHIRLY1 contributes significantly to pathogen resistance in SA-dependent lines [110].
The only one identified NF-Y was down-regulated (Fig 7). It has shown that an increasing number of NF-Ys play significant roles in the various stages of root node coexistence and arbuscular mycorrhizae as well as in the interaction of plants with pathogenic microorganisms [111]. Barley NF-Ys comprise a relatively diverse gene family involved in various processes such as starch metabolism, sucrose metabolism, and photosynthesis. There is less evidence about the role of NF-Ys against pathogenic attacks. Overexpression of rice NF-YA (OsHAP2E) causes resistance to Magnaporthe oryzae and Xanthomonas oryzae [112]. The signaling pathways of the NF-YAs regulate beneficially host-microorganism interactions and indicate the existence of common mechanisms between both types of biological interactions. The results demonstrate that some members in the miR169 family might handle NF-Ys to activate the ISR pathway [113]. Using these findings, we can illustrate the contrast between regulators allowing plants to construct complex networks of interacting pathways in the presence of various pathogens and hormonal signals.
The GNAT family is generally considered to be comprised of four subfamilies designated GCN5, ELP3, HAT1, and HPA2. It has shown that members of the GNAT superfamily activate gene transcription by the promotion of euchromatin formation and exert various levels of contributions to cell morphogenesis and virulence [114]. In addition, a GNAT acyltransferase has been identified that plays a role in multiple biological processes, including the regulation of growth, morphogenesis, cellulase expression, and pathways related to biotic stress responses and hormone signaling. For example, in Beauveria bassiana, three HATs, including a GNAT superfamily member, have previously been characterized to variously affect developmental and virulence pathways in this fungus.
In total, analyses revealed that most TF families such as AP2/ERF, bHLH, bZIP, WHIRLY1, and NF-Ys are not only implicated in the regulation network of SAR or ISR but also are involved in the regulation of SA and ET/JA cross-talks [115]. Transcriptional regulator NPR1 has been reported to be required for basal resistance toward biotrophic and hemibiotrophic pathogens, as well as the development of SAR and ISR [37]. SA-dependent SAR requires a nuclear function of NPR1, while SA-independent ISR requires cytosolic [116,117] NPR1 localization and involves the JA and ET [118]. Therefore, post-translational modifications of NPR1 appear to play a significant role in modifying transcriptional defense outputs. In particular, we found that WHIRLY1 and NF-Ys in interaction with the NPR1 regulator play a more important role in connection with SAR/ISR pathways, respectively, and hormone-induced defenses thus can be the center of attention and more study.

PKs in regulation of signaling networks and immune response
PKs are the largest family of molecular switches, which can regulate protein activities and basic cellular functions. However, only a fraction of PKs is functionally characterized even in model plant species. To investigate signal transduction mechanisms in this study, we identified potential PKs that play key roles in signaling networks and plant responses to various biotic stresses and hormonal signaling. We detected 3 PKs among the biotic DEGs and 3 PKs among the hormonal DEGs that are classified into RLK/Pelles and TKL groups (Table 4). RLK/Pelles play important roles from growth regulation to defense response [119]. The dramatic expansion of this family has been considered very important for plant-specific adaptations [119]. We realized that hundreds of RLK/Pelles are up-regulated by biotic stresses. In addition, the rate of stress responsiveness is related to the degree of consecutive duplication in RLK/Pelle subfamilies [120]. Members of this family are involved in many different plant processes, including regulation of meristem proliferation, reproduction, hormonal signal transduction, and response to biotic and abiotic stresses [121,122]. Accordingly, these results suggest that PKs play a central role in signal transduction during plant growth and development and also in plant responses to biotic stresses. One-third of the barley RLK/Pelle members were down-regulated (Table 4) during development, suggesting that these genes may have a negative effect on regulatory functions. Many non-RLK groups, including all of TKLs, however, exhibited high expression levels during development, confirming that these genes may play positive regulatory roles in plant growth. This result shows a positive or negative regulatory function among PKs.
Therefore, RLK/Pelles and TKLs evolved as the families of highly dynamic molecular switches and components of signaling networks, which leads to constitutive activation of the defense response, accumulation of SA, and increase of NPR1 performance [123]. The accumulation of SA can further induce SAR by up-regulating PRs. In particular, the fact that the most barley RLK families were up-regulated by SA, which are involved in SAR-related gene regulation, makes them interesting candidates for bioengineering broad-spectrum resistance genotypes.

miRNAs in regulation of gene expression
Plants have evolved complex mechanisms involving a diversity of small RNA-guided stress regulatory networks that can provide new insights for genetically improved plant pathogen resistance. MicroRNAs (miRNAs) are a major class of small non-coding RNAs that have emerged as significant regulators of gene expression during or after transcription [124]. Some miRNAs have been related to pathogen responses, and they play important roles in plants infected by bacteria, viruses, nematodes, and fungi [125].
Barley is one of the most important crops that has been considered for miRNA studies [124]. In the present study, 432 biotic miRNAs belonging to 56 families and 85 hormonal miR-NAs belonging to 32 families, and also 15 common miRNAs for both biotic and hormonal data were identified (S6 Table). A high proportion of these miRNAs are involved in the different BPs (Fig 8B). The major members of the miRNAs for both studies belong to conserved families, including miR6192, miR6180, and miR6214, miR6196, miR6184, miR166, and miR156. Based on a research by Djami-Tchatchou et al. (2017) about 100 miRNAs were identified for the first time in barley using deep sequencing [124]. As compared these miRNAs to our results, we confirmed that many microRNAs (miR156, miR159, etc.) showed orthologs in wheat or rice, while it appears that several microRNAs (miR444) were expressed in barley specifically.
A combination of small RNA and mRNA degradome analyses to recognition and characterization of many conserved miRNAs (miR156, miR166, etc.) and novel miRNAs (hvu-miR1120b) together with several miRNA target genes were done (Fig 8B). The results showed that some miRNAs are involved in different functions, including carbohydrate translocation, cell differentiation, defense response, photosynthesis, and phytohormone signaling pathways [126]. This finding revealed that miRNAs play a significant role in regulating early growth, the development of barley, and the signaling pathways of different defense hormones [127].
Since ISR differs from SAR in stimulating inducers and signal transduction, different miR-NAs may participate in ISR and SAR. miR393 was the first identified miRNA responsible for triggering the SAR and opened a way to study the functions of miRNAs in SAR [113]. Some miRNAs like miR472, miR482, and miR2109 suppress resistance (R) genes to regulate defense responses [128]. Although miR846 and miR825 have been reported to regulate Bacillusinduced ISR, they belong to non-conserved miRNAs, which exist only in some plant species [129]. Bacillus velezensis FZB42 can activate ISR to enhance plant defense response against pathogen infections. Two miRNAs, miR169 and miR395, which were repressed in FZB42-treated plants especially, were selected as candidate ISR-associated miRNAs. All of them contained at least one cis-element associated with defense response, whereas the miR169 family had the most abundant defense-related cis-elements [113]. In addition to the regulation of some abiotic stresses, the miR169 family and its targets NF-YA genes participate in plant immunity [111]. This miRNA acts as a negative regulator in rice immunity against the blast fungus Magnaporthe oryzae by repressing the NF-YAs [130].
In summary, our results showed that some members in the miR169 family ( Fig 8A) might regulate NF-Y TFs to activate the ISR. This study is helpful in better understanding the regulatory roles of barley miRNAs in hormonal signaling and pathogen-activated ISR.

Promoter analysis and potential regulatory elements associated with the pathogen resistance pathways
Promoters are regions of the genome that can bind TFs. The study of promoters, TFs, and downstream genes is a fascinating topic in post-genomics and can provide new insights into vital BPs [131,132]. We performed this analysis to identify cis-acting elements located upstream of DEGs according to whether these are under common regulatory control. Eleven motifs with significant scores were identified (S8 Table).
Many cis-regulatory elements are responsive to pathogens and hormones. The most frequent cis-acting elements were highly matched to the MA1267.1 motif (Tables 5 and 6 and S7  Table), a C2H2 zinc finger motif (DOF5.8). These TFs play an important role in many BPs, e.g., flowering time, seed development, and responses to hormones and stresses [133]. He et al. (2015) [134] realized that the Arabidopsis DOF5.8 is an upstream regulator of a gene encoding an NAC family member in response to biotic and abiotic stresses [135]. More particularly, DOF5.8 was shown to play a significant role in the correct development of veins within the leaves of A. thaliana, and its promoter was demonstrated to be regulated by auxin [136]. Actually, the overexpression of this TF leading to a modification within the expression of the many genes, was involved in auxin biosynthesis, plant tissue formation, secondary cell wall deposition and similarly as signaling during physiological processes [137,138].
Since cross-talk between plant hormone signaling pathways is critical for controlling the immune response during pathogen invasion, we investigate a basic helix-loop-helix (bHLH) transcription activator, which activates JA signaling in cereal when is localized to the nucleus [139]. The bHLHs are one of the largest TF families and participate in the regulation of plant growth and development [140,141], iron homeostasis [142], root vascular cell proliferation [143], shoot branching [144], stomatal initiation [145], flowering time [146], pollen, gynoecium and fruit development [147], and grain yield [148]. bHLH is involved in the phytohormone-mediated pathogen stress adaptation and resistance development [149]. Following infection by fungal pathogens such as Magnaporthe oryzae, the SA signaling regulator, NPR1 sequesters bHLH in the cytosol, which activates SA signaling but represses JA signaling. We, therefore, propose that bHLHs may regulate the SA/JA antagonism in cereals [139]. Further, wheat bHLH overexpression negatively regulated resistance to Pseudomonas syringae through JA and ET signaling in transgenic Arabidopsis [150]. Although some functions of bHLHs have been characterized, the biological functions of most of them remain unclear, especially in gramineous crops such as rice, maize, and wheat [149]. The subgroup of bHLHs (bHLH3, bHLH13, bHLH14, and bHLH17) in Arabidopsis was found to act as targets of JAZs and negatively regulate JA-mediated defense and development [151]. bHLHs are repressors that bind to cis-regulatory elements of MYC target genes and inhibit transcription [152]. In other words, stabilized JAZs and bHLH repressors counteract the activation of JA responses to fine-tune defense and development. Therefore, the accumulation of SA in response to biotrophs and insects interferes with JA defense gene expression via other routes [152].

Co-expression network analysis and identification of hub genes
Co-expression network analysis of the DEGs was performed to discover processes and functions involved in response to pathogens and hormones (Fig 5 and S3 Table). Interestingly, 4 modules, including brown, blue, and turquoise for both datasets, were significantly related with response to defense pathways. After studying the genes within these modules, finally, four genes (SSI2, PR1, PKT3, and LOX2) related to biotic studies and two genes (AOS and OPR3) related to hormonal studies that were involved in the regulation of BPs were identified to further investigation (Fig 9 and S10 Table).
The Arabidopsis ssi2 is a mutant line affected in both SA-and JA-mediated defense signaling. SSI2 encodes a plastid-localized stearoyl-acyl carrier protein desaturase (SACPD) that desaturates stearic acid to oleic acid in the chloroplast. The ssi2 line causes a truncation in the SSI2, which results in the loss of 90% of SACPD activity [153]. SACPD is an archetypical member of a family of soluble FAs desaturases; these enzymes play an important role in regulating the overall level of desaturated FAs in the cell and in modulating cross-talk among different defense signaling pathways [154]. The ssi2 lines are stunted in size, show spontaneous death of the cell, constitutive overexpress of PR1, accumulation of high levels of SA, and exhibit enhanced resistance to pathogens. Taken together, SSI2 appears to affect the activation of several defense responses by modulating a NIM1/NPR1-independent defense pathway. The double mutant lines of ssi1 npr1-5 and ssi2 npr1-5, suppressed in constitutive SA signaling, show necrotic lesions and basal-level expression of PR1 in the absence of a functional NIM1/NPR1. Therefore, an NPR1-dependent signaling pathway acts additively with the ssi2-induced NPR1-independent signaling pathway to enhance PR1 expression [155]. In addition to regulating SAR, the NIM1/NPR1 plays a role in ISR as well [156], which is activated in Arabidopsis after exposure of roots to certain nonpathogenic Pseudomonas fluorescens strains. Unlike SAR, induction of ISR does not require SA accumulation and is not related to the expression of SAR-associated PRs. However, like SAR, induction of ISR is dependent upon NIM1/NPR1 [156]. Therefore, SAR and ISR appear to represent distinct defense pathways in Arabidopsis that are likely to differ in their modes of action, yet share a requirement for NIM1/NPR1 [157]. How NIM1/NPR1 modulates the expression of both SAR and ISR is unclear, but it is becoming apparent that NIM1 plays a central role in multiple defense signaling pathways [158]. Because SACPD catalyzes a desaturation step required for JA biosynthesis [153], we attempted to determine whether the treatment of barley seedlings by MeJA can affect the SSI2 expression ( Fig 11A).
PR proteins, which are classified into 17 families, accumulate after pathogen infection or related situations in many plant species and cause defense responses, especially during the hypersensitive response (HR) [159]. Among the PR family, PR1 has been frequently applied as a marker gene for SAR in many species [160]. One PR1 homolog, At2g14610, encodes a basic protein that has been identified in many plant species as a marker of defense against oomycetes, fungi, bacteria, or viruses [161]. From these results, we tend to suppose that only PR1 relates directly to pathogen resistance in Arabidopsis and the others contribute to other functions. In response to defense hormones, PR1 was induced briefly by MeJA treatments ( Fig  11B). Indeed, the MeJA treatment and pathogen infection induce the same set of PRs, suggesting JA/ET-signaling pathways are involved in pathogen resistance [162,163].
Another hub gene is PKT3 which encodes a 3-ketoacyl thiolase in peroxisomes that catalyzes the last step of the β-oxidation of fatty acids to produce one molecule of acetyl-CoA in each repeat cycle. This gene shows higher expression than the redundant homologs, such as PKT1 and PKT5. It is now broadly accepted that PKT3 acts in IAA biosynthesis through its βoxidative decarboxylation in peroxisomes [164]. The coordinated induced expression of β-oxidation genes is essential to provide the energy supply for germination or post-development of germination. Mechanical damage triggers the local and systemic induction of PKT3 in Arabidopsis. Although most of the β-oxidation genes were activated by wound-related factors such as dehydration and ABA, JA induces only ACX1 and PKT5 [101]. Reduced expression of ACX1 or PKT3 in transgenic plants, expressing their corresponding mRNAs in antisense orientation, was correlated with defective wound-activated synthesis of JA and reduced expression of JA-responsive genes. Induced expression of JA-responsive genes by exogenous application of JA was unaffected in those transgenic plants, suggesting that ACX1 and PKT3 play a major role in driving wound-activated responses by participating in the biosynthesis of JA [101].
Another hub gene is AOS, a key enzyme of the JA signaling pathway. The mutation of the single AOS in Arabidopsis leads to the complete elimination of JA production [81]. It has proved that the plants constitutively expressing AOS delay pathogen progression and exhibit more vigorous growth. They also showed the up-regulation of JA and SA-related genes, increased resistance to pathogens, and robust induction of PRs, such as PR1, PR3, and PR5 [165,166]. Taken together, these results demonstrate that the AOS is an important player in plant defense responses against pathogen infections and affects the JA biosynthetic pathway thus its expression in susceptible varieties may be a valuable tool to mitigate biotic stress responses [167].
Arabidopsis OPDA (12-oxophytodienoic acid) reductase (OPR3) catalyzes the reduction of OPDA to 3-oxo-2-(2 0 (Z)-pentenyl)-cyclopentane-1 octanoic acid (OPC-8:0). OPDA is known as a biosynthetic precursor of the JA. Although Arabidopsis contains at least two other OPR homologs, OPR1 and OPR2, their protein products do not catalyze the reduction of OPDA [168]. Among the 6 wound-induced OPRs of A. thaliana, only OPR3 is involved in JA biosynthesis [169]. Interestingly, inducing OPR3 expression by JA indicates an opportunity for feedback regulation of gene expression. Expression of opr3 in plants treated with OPDA differs from that in opr3 lines treated with JA. In contrast to OPDA, JA induces the transcription of three homologous genes. Research results show that OPDA has different regulatory functions than JA and is an independent signaling molecule [170]. Generally, the plants are deficient in the biosynthesis of JA, and they accumulate OPDA when wounded. Significantly, OPR3 has capable JA defense responses against insect pests and pathogens such as Botrytis cinerea and Alternaria brassicicola. Although topical application of JA and MeJA induces transcription of defense genes in Arabidopsis, OPDA alone is sufficient for these responses. This raises the question of which of these hormones represent the active signal molecules in plants [169]. The OPR3 is the candidate hub gene that is validated for qPCR.
Lipoxygenase (LOX) is last hub gene candidate for more validation in qPCR. An enzyme that catalyzes the hydro-peroxidation of specific unsaturated fatty acids as well as initiates JA synthesis and converts linoleic acid into 13-hydroperoxy linolenic acid Therefore, LOXs are induced in the ISR pathway and involved in defense responses against the pathogen [171]. Among the 6 LOXs of Arabidopsis, four of them are 13-LOXs (LOX2, LOX3, LOX4, LOX6), although their functions are only partly understood. It was thought that LOX2 was involved in wound response for a long time [169]. The subsequent studies revealed that LOX2 was responsible for the bulk of JA formation in pollen and stamen development in the first hours of wounding [172]. In mutants with defects in a JA-regulated signal transduction pathway, the regulation of JA-inducible genes might be altered. It has been reported that A. thaliana LOX2 is induced by MeJA. LOX2 is also involved in lipid peroxidation that occurs under abiotic and biotic stresses. A LOX2 mediated double oxygenation of plastid galactolipids leading to the production of arabidopsides was observed upon pathogen infection, however it was not responsible for the pathogen-induced increase in JA [173]. Interestingly, the formation of lipid peroxides accompanies the synthesis of azelaic acid, a new signaling compound priming the immune responses [174].

Validation of candidate genes
The results suggest that the gene expressions are reliable and mainly are consistent with those from data analysis. After MeJA treatment, as shown in Fig 11, the expression of SSI2 was early induced. This gene causes biotic stress-induced constitutive SA signaling [153]. High levels of linolenic acid in ssi2 mutant lines are required for the activation of certain ET/JA-mediated responses, because these lines are responsive to both of these hormones, but show normal levels of SA and PRs expression [153]. JA treatment would not be sufficient to activate PDF1.2 expression in ssi2 lines or restore resistance to pathogens since these mutants lack or have low levels of this co-activating signal [154]. Supporting this possibility is the discovery that the induction of PR1 or injection of oleic acid into the leaves of these lines restores JA-inducible PDF1.2 expression [175].
The gene expression analysis indicated that PKT3 and PR1 were induced to a significantly higher level in response to MeJA treatment (Figs 11 and S6). The induced PR1 is a useful molecular marker for the SAR responses and is essential for the transduction of the SA signals [176]. The induction of the same set of PRs in response to pathogen infections, MeJA or ethephon, suggests the involvement of JA/ET-signaling pathways in mediating resistance against pathogens consistent with the potential hemibiotrophic nature. Indeed, the responses of the PR1 homologs to defense signaling compounds indicate that the signaling pathways in rice conferred by SA and JA were synergistic as well as that by ET and JA, and SA and ET [160].
The induction of LOX2, OPR3, and AOS activates the defense responses against the pathogens and biosynthesis of JA, thus they can induce the ISR pathway [81]. These genes involved in JA biosynthesis were significantly induced during Botrytis cinerea infection in wild type, but their expression was reduced in the coi1 mutant that is a critical component of the JA receptor [177]. Additionally, the activation of these genes was also repressed by SA, suggesting that the JA biosynthesis pathway may be a target for SA-mediated antagonism [178]. According to these results, based on global differential expression matrices, machine-learning approaches could help predict how potential plant elicitors act.

Conclusions
Using barley transcriptome data, we integrated a meta-analysis method with gene co-expression network analysis. Based on these analyses, we identified genes that play a significant role against pathogens and in hormonal signaling. The results indicate that most of the DEGs influence a wide range of biological processes related to pathogen resistance and hormone signaling transduction. Because of the challenging environment, crops encountered numerous biotrophic and necrotrophic pathogens that required SA and JA signals to respond. The importance of TFs (e.g., WRKYs, MYBs, bHLHs, and bZIPs), kinases (e.g., RLK/Pelles), and miRNAs (e.g., miR6192) in the responses to pathogens was demonstrated and several novel interesting genes with unknown functions were identified. We showed that WHIRLY1 and NF-Y TF families, and miR169 family have a more important role in connection with SAR /ISR pathways and hormone-induced defenses. The results can be the center of more attention and study. The promoter analysis identified 11 overrepresented cis-acting elements such as zinc finger proteins in upstream regions of DEGs, which improve plant defense against fungal attack by ISR-dependent pathway. Moreover, it confirmed the results of the meta-analysis by showing the existence of pathogen-specific and consensus modules against pathogens and in hormonal signaling. The PPI network analysis helped to visualize and introduce the key genes such as SSI2, PR1, PKT3, AOS, LOX2, and OPR3 related to hormonal signaling and pathogenic infections that were validated in qPCR. According to integrative data analysis, these genes showed a significant expression change after MeJA treatment in qPCR, the potential candidates for the regulation of JA biosynthesis. A number of the DEGs with high connectivity and conserved expression but with weak annotation were also recognized. It may contribute to a much better understanding of how barley responds to MeJA treatment resulting in ISR. Unlike SAR, induction of ISR does not require SA accumulation and is not related to the expression of SAR-associated PRs. However, like SAR, induction of ISR is dependent upon NIM1/NPR1. Thus, in addition to regulating SAR, NPR1/NIM1 has also been shown to be involved in the activation of ISR by the SSI2. We suggest these genes for more study, which might provide a more robust bio-signature for phenotypic traits as a more promising resource of molecular biomarkers and potential candidate genes for engineering pathogen resistance in barley. The different letters on the columns represent the significant differences given by Duncan's multiple range statistical analysis. Vertical bars indicate ± SE of the mean (n = 3, P-value � 0.05). The X-axis shows the times after MeJA treatment. As it is known, the target genes showed the highest expression at the time points of 6 h to 24 h and showed the lowest expression at the time of 48 h. (TIF) S1